<<<<<<< HEAD

Εργαστηριακή Άσκηση 2

Σκοπός της δεύτερης σειράς ασκήσεων είναι η εξοικείωση με τις συναρτήσεις σχεδιασμού φίλτρων πεπερασμένης κρουστικής απόκρισης (FIR) και την υλοποίησή τους στο MATLAB. Προτού ξεκινήσετε την άσκηση θα πρέπει να μελετήσετε με προσοχή το Κεφάλαιο 1 και, ειδικότερα, την παράγραφο 1.3 του τεύχους του μαθήματος. Το MATLAB (www.mathworks.com) είναι ένα διαδραστικό εμπορικό πρόγραμμα (Windows, Linux, Unix) με το οποίο μπορείτε να κάνετε εύκολα αριθμητικές πράξεις με πίνακες. Μπορεί να το έχετε εγκατεστημένο τοπικά, στον προσωπικό σας υπολογιστή ή να εργάζεστε σε κάποιο Εργαστήριο Προσωπικών Υπολογιστών (ΕΠΥ) της Σχολής σας που διαθέτει το συγκεκριμένο λογισμικό.

Live Code

Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).

Python Interactive Code

Μέρος 1: Εισαγωγή

Στο MATLAB οι συναρτήσεις \(fft\) και \(ifft\) υποθέτουν ζεύγος μετασχηματισμού Fourier \(x(t)\) και \(X(f)\) υπολογισμένων σε μη αρνητικά διαστήματα \( t=[0:N-1]ts \) και \( f=[0:N-1]fo. \) Όπως έχετε ήδη δει στην Εργαστηριακή άσκηση 1, το άνω μισό μέρος του διαστήματος συχνοτήτων αντιστοιχεί στις αρνητικές συχνότητες του σήματος, όταν υπολογίζουμε το \(X(f)\) με τη βοήθεια της συνάρτησης \(fft\) (για σήματα πραγματικών τιμών, αυτά τα δυο μισά είναι κατοπτρικά ως προς το μέσο του διαστήματος). Ακριβώς το ίδιο ισχύει και για το άνω μισό μέρος του χρονικού διαστήματος, όταν το σήμα \(x(t)\) προκύπτει από τον αντίστροφο μετασχηματισμό Fourier μέσω της \(ifft\).

Το MATLAB διαθέτει τη συνάρτηση \(fftshift\) για να ολισθήσει κυκλικά τις τιμές του σήματος ή του μετασχηματισμού Fourier, ώστε να αντιστοιχούν σε κεντραρισμένα στο μηδέν αμφίπλευρα διαστήματα, δηλαδή, στις χρονικές στιγμές \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) ή στις συχνότητες \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). Με τον τρόπο αυτό μπορούμε να παράγουμε τα \(xb(t)\) και \(Xb(f)\) που αντιστοιχούν στην αμφίπλευρη αναπαράσταση του σήματος και του μετασχηματισμού Fourier.

Για να κατανοήσετε τα ανωτέρω θεωρείστε το διάνυσμα \([1 2 3 4]\) ως το αποτέλεσμα του \(FFT\) μήκους 4. Τότε, το πρώτο στοιχείο (1) είναι ο όρος dc, το τρίτο στοιχείο (3) είναι το σημείο στο μισό της συχνότητας δειγματοληψίας \(fs/2\), που μπορεί να εκληφθεί ότι αντιστοιχεί είτε στην \(–fs/2\) είτε στην \(+fs/2\). Τα στοιχεία 2 και 4 αντιστοιχούν στις συχνότητες \(+fs/4\) και \(–fs/4\). Εφαρμόζοντας την \(fftshift\), το στοιχείο 3 εμφανίζεται πρώτο, που σημαίνει ότι στο MATLAB αντιστοιχεί στην αρνητική συχνότητα \(–fs/2\), το επόμενο στοιχείο 4 αντιστοιχεί στη συχνότητα \(–fs/4\) ακολουθούμενο από το dc και τη συχνότητα \(+fs/4\). Για ένα μετασχηματισμό περιττού μήκους, δεν υφίσταται σημείο για το \(±fs/2\). Έτσι για το διάνυσμα \([1 2 3]\), η εφαρμογή της \(fftshift\) θα δώσει τα στοιχεία που αντιστοιχούν στις συχνότητες \(–fs/3, 0, +fs/3\).

Εκτός του ότι παράγουν εξόδους με τις αρνητικές συχνότητες ή χρόνους στο άνω μισό του διανύσματος, αμφότερες οι συναρτήσεις \(fft\) και \(ifft\) αναμένουν ως είσοδο διάνυσμα με την ίδια μορφή, αφού προφανώς ισχύουν οι ταυτότητες

\(h == \text{ifft}(\text{fft}(h))\) και \(H == \text{fft}(\text{ifft}(H))\)

Η πρώτη υποδεικνύει ότι η είσοδος της \(ifft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(fft\), και η δεύτερη ότι η είσοδος της \(fft\) πρέπει να είναι αντεστραμμένη, όπως την παράγει η \(ifft\).

Στο επόμενο σχήμα βλέπετε παραστατικά ένα ημιτονικό σήμα που έχει πολλαπλασιασθεί με παράθυρο Blackman τόσο στην αμφίπλευρη, όσο και την μονόπλευρη αναπαράστασή του.

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import Layout
from IPython.display import display, clear_output

np.random.seed(0)

# Generate random data
t = np.linspace(-40, 40, 80)  # 80 samples
random_data = np.random.randn(len(t))

# Initial start and end points for the window
start, end = 30, 50

# Output widget for the plots
output = widgets.Output()

def apply_window_and_plot(start, end):
    with output:
        clear_output(wait=True)
        # Create a sinusoidal window
        sin_window = np.zeros_like(t)
        sin_window[start:end] = np.sin(np.linspace(0, np.pi, end - start))

        # Apply sinusoidal window to random data
        windowed_data = sin_window * random_data

        # Append the negative times to the end of the time array
        t_append_negative = np.concatenate((t[t >= 0], t[t < 0] + t[-1] - t[0] + (t[1] - t[0])))
        windowed_data_append_negative = np.concatenate((windowed_data[t >= 0], windowed_data[t < 0]))
        plt.close('all')  # Close all existing figures
        # Create subplots
        fig, axs = plt.subplots(2, 1, figsize=(10, 10))

        # Plot original windowed data
        axs[0].stem(t, windowed_data, linefmt='C0-', markerfmt='C0o', basefmt=" ")
        axs[0].set_title("Original Sinusoidal Windowed Random Data", fontsize=16)
        axs[0].set_xlabel("Time (samples)", fontsize=12)
        axs[0].set_ylabel("Amplitude", fontsize=12)

        # Plot with negative times appended
        axs[1].stem(t_append_negative, windowed_data_append_negative, linefmt='C1-', markerfmt='C1o', basefmt=" ")
        axs[1].set_title("Sinusoidal Windowed Random Data with Negative Times Appended", fontsize=16)
        axs[1].set_xlabel("Time (samples)", fontsize=12)
        axs[1].set_ylabel("Amplitude", fontsize=12)

        # Adjust layout and show plot
        plt.tight_layout()
        
        plt.show()
        

# Define the range slider
window_range_slider = widgets.IntRangeSlider(
    value=[start, end],
    min=0,
    max=len(t),
    step=1,
    description='Window Range:',
    layout=Layout(width='90%'),  # Setting the width to 90%
    continuous_update=False
)

# Function to update the plot based on the range slider
def on_range_change(change):
    new_start, new_end = change['new']
    apply_window_and_plot(new_start, new_end)

# Observe the range slider for changes
window_range_slider.observe(on_range_change, names='value')

# Group the slider and output widget in a vertical box layout
vbox = widgets.VBox([window_range_slider, output])

# Display the VBox
display(vbox)

# Initialize with the default plot
apply_window_and_plot(start, end)

<<<<<<< HEAD

=======

>>>>>>> d546939e94d695c6e9c44c085bc5ec2e2e14440e =======

Lab Exercise 2

The purpose of the second series of exercises is to familiarize with the design functions of finite impulse response (FIR) filters and their implementation in Python & MATLAB.

Tip

Before you start the exercise, you should carefully study Chapter 1 and, in particular, paragraph 1.3 of the course booklet.

>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Part 1: Training

Live Code

Press the following button to make python code interactive. It will connect you to a kernel once it says “ready” (might take a bit, especially the first time it runs).

Python Interactive Code

Part 1: Introduction

In MATLAB, the functions \(fft\) and \(ifft\) assume a pair of Fourier transforms \(x(t)\) and \(X(f)\) calculated in non-negative intervals \(t=[0:N-1]ts\) and \(f=[0:N-1]fo. \) As you have already seen in Laboratory Exercise 1, the upper half of the frequency interval corresponds to the negative frequencies of the signal, when calculating \(X(f)\) with the help of the \(fft\) function (for real-valued signals, these two halves are mirror images with respect to the middle of the interval). The same applies to the upper half of the time interval, when the signal \(x(t)\) results from the inverse Fourier transform through \(ifft\).

MATLAB has the \(fftshift\) function to cyclically shift the values of the signal or the Fourier transform so that they correspond to zero-centered bilateral intervals, i.e., at times \(tb=[-ceil((N-1)/2): floor((N-1)/2)]ts\) or frequencies \(fb=[–ceil((N-1)/2): floor((N-1)/2)]fo\). This way, we can produce \(xb(t)\) and \(Xb(f)\) that correspond to the bilateral representation of the signal and the Fourier transform.

To understand the above, consider the vector \([1 2 3 4]\) as the result of a length 4 \(FFT\). Then, the first element (1) is the dc term, the third element (3) is the point at half the sampling frequency \(fs/2\), which can be interpreted as corresponding either to \(–fs/2\) or \(+fs/2\). The elements 2 and 4 correspond to frequencies \(+fs/4\) and \(–fs/4\). By applying \(fftshift\), the element 3 appears first, meaning that in MATLAB it corresponds to the negative frequency \(–fs/2\), the next element 4 corresponds to the frequency \(–fs/4\) followed by the dc and the frequency \(+fs/4\). For an odd-length transformation, there is no point for \(±fs/2\). Thus, for the vector \([1 2 3]\), applying \(fftshift\) will yield the elements corresponding to the frequencies \(–fs/3, 0, +fs/3\).

Besides producing outputs with negative frequencies or times in the upper half of the vector, both \(fft\) and \(ifft\) functions expect as input a vector in the same format, since obviously the identities hold

\(h == \text{ifft}(\text{fft}(h))\) and \(H == \text{fft}(\text{ifft}(H))\)

The first indicates that the input of \(ifft\) should be reversed, as produced by \(fft\), and the second that the input of \(fft\) should be reversed, as produced by \(ifft\).

In the next figure, you see illustratively a sinusoidal signal that has been multiplied by a Blackman window both in the bilateral, and the unilateral representation.

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from ipywidgets import Layout
from IPython.display import display, clear_output

np.random.seed(0)

# Generate random data
t = np.linspace(-40, 40, 80)  # 80 samples
random_data = np.random.randn(len(t))

# Initial start and end points for the window
start, end = 30, 50

# Output widget for the plots
output = widgets.Output()

def apply_window_and_plot(start, end):
    with output:
        clear_output(wait=True)
        # Create a sinusoidal window
        sin_window = np.zeros_like(t)
        sin_window[start:end] = np.sin(np.linspace(0, np.pi, end - start))

        # Apply sinusoidal window to random data
        windowed_data = sin_window * random_data

        # Append the negative times to the end of the time array
        t_append_negative = np.concatenate((t[t >= 0], t[t < 0] + t[-1] - t[0] + (t[1] - t[0])))
        windowed_data_append_negative = np.concatenate((windowed_data[t >= 0], windowed_data[t < 0]))
        plt.close('all')  # Close all existing figures
        # Create subplots
        fig, axs = plt.subplots(2, 1, figsize=(10, 10))

        # Plot original windowed data
        axs[0].stem(t, windowed_data, linefmt='C0-', markerfmt='C0o', basefmt=" ")
        axs[0].set_title("Original Sinusoidal Windowed Random Data", fontsize=16)
        axs[0].set_xlabel("Time (samples)", fontsize=12)
        axs[0].set_ylabel("Amplitude", fontsize=12)

        # Plot with negative times appended
        axs[1].stem(t_append_negative, windowed_data_append_negative, linefmt='C1-', markerfmt='C1o', basefmt=" ")
        axs[1].set_title("Sinusoidal Windowed Random Data with Negative Times Appended", fontsize=16)
        axs[1].set_xlabel("Time (samples)", fontsize=12)
        axs[1].set_ylabel("Amplitude", fontsize=12)

        # Adjust layout and show plot
        plt.tight_layout()
        
        plt.show()
        

# Define the range slider
window_range_slider = widgets.IntRangeSlider(
    value=[start, end],
    min=0,
    max=len(t),
    step=1,
    description='Window Range:',
    layout=Layout(width='90%'),  # Setting the width to 90%
    continuous_update=False
)

# Function to update the plot based on the range slider
def on_range_change(change):
    new_start, new_end = change['new']
    apply_window_and_plot(new_start, new_end)

# Observe the range slider for changes
window_range_slider.observe(on_range_change, names='value')

# Group the slider and output widget in a vertical box layout
vbox = widgets.VBox([window_range_slider, output])

# Display the VBox
display(vbox)

# Initialize with the default plot
apply_window_and_plot(start, end)

The corresponding Fourier transforms, in bilateral and unilateral representation, are shown in the next figure:

When \(x(t)\) and \(X(f)\) are produced by MATLAB, no special attention is needed, except for the cyclic shift in case, for example, we want to plot the bilateral spectrum or signal. However, when one of \(x(t)\) or \(X(f)\) is defined by the user, more attention is required, because usually, the bilateral signals or spectra are used. You can move from one representation to the other as follows:

\(x = ifftshift(xb)\), \(X = fft(x)\), \(Xb = fftshift(X)\), if you start from a bilateral signal and want to end up with a bilateral spectrum, and

\(X = ifftshift(Xb)\), \(x = ifft(X)\), \(xb = fftshift(x)\), if you start from a bilateral spectrum and want to end up with a bilateral signal,

where the MATLAB function \(ifftshift\) performs the inverse operation of \(fftshift\). When \(N\) is even, \(fftshift\) and \(ifftshift\) yield the same result. However, when \(N\) is odd, this is not the case, and care must be taken in their use. In practice, careful application of the above is important when calculating the phase of the spectrum. The amplitude of the spectrum is not affected by the cyclic shift of elements caused by \(fftshift\) and \(ifftshift\) (see DFT properties).

Exercise

Try the following in the command window to consolidate the use of the \(fftshift\) and \(ifftshift\) functions.

Python

import numpy as np

# Create the array X
X = np.arange(-2, 3)  # equivalent to MATLAB's -2:2

# Apply fftshift and ifftshift
fftshifted_X = np.fft.fftshift(X)
ifftshifted_X = np.fft.ifftshift(X)

# Double fftshift and a combination of fftshift and ifftshift
Y = np.fft.fftshift(np.fft.fftshift(X))
Z = np.fft.ifftshift(np.fft.fftshift(X))

# Check if arrays are equal
X_equals_Y = np.array_equal(X, Y)
X_equals_Z = np.array_equal(X, Z)

print(f"X equals Y: {X_equals_Y}")
print(f"X equals Z: {X_equals_Z}")

Matlab

X = [-2:2]
fftshift(X)
ifftshift(X)
Y = fftshift(fftshift(X));
Z = ifftshift(fftshift(X));
isequal(X,Y)
isequal(X,Z)

X equals Y: False
X equals Z: True

<<<<<<< HEAD

Ερώτηση 1: Ποιο εκ των διανυσμάτων Υ και Ζ ισούται με το X; Γράψτε την απάντησή σας σε ένα αρχείο κειμένου lab2_nnnnn.txt, όπου nnnnn τα πέντε τελευταία νούμερα του αριθμού μητρώου σας, χρησιμοποιώντας το Notepad από το μενού των Windows (Start → Programs → Accessories → Notepad) και αποθηκεύστε το στον φάκελο My Documents. Θα υποβάλετε το αρχείο αυτό ηλεκτρονικά στο τέλος, αφού απαντήσετε και τις επόμενες ερωτήσεις, οπότε μπορείτε να τα αφήσετε ανοικτό.

Ερώτηση 2: Επαναλάβατε με \(X=[-1:2]\). Τι παρατηρείτε; Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt. Δοκιμάστε στο παράθυρο εντολών τα ακόλουθα δύο παραδείγματα για να εμπεδώσετε τη χρήση των συναρτήσεων \(fftshift\) και \(ifftshift\) σε συνδυασμό με τις \(fft\) και \(ifft\).

=======

Question 1: Which of the vectors Y and Z equals X? Write your answer in a text file lab2_nnnnn.txt, where nnnnn are the last five numbers of your registration number, using Notepad from the Windows menu (Start → Programs → Accessories → Notepad) and save it in the My Documents folder. You will submit this file electronically at the end, after answering the next questions, so you can leave it open.

Question 2: Repeat with \(X=[-1:2]\). What do you observe? Write your answer in the text file lab2_nnnnn.txt. Try the following two examples in the command window to consolidate the use of the \(fftshift\) and \(ifftshift\) functions in combination with \(fft\) and \(ifft\).

>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Code

Python

import numpy as np
from scipy.fft import fft, ifft, fftshift, ifftshift
import matplotlib.pyplot as plt
from ipywidgets import interact, Layout
import ipywidgets as widgets

# First part: Original signal and its FFT
xb = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1])
x = ifftshift(xb)
X = fft(x)
Xb = fftshift(X)  # Spectrum with DC component in the center

# Second part: Low-pass filter effect
Xb_low_pass = np.array([0, 0, 1, 1, 1, 1, 1, 0, 0])
X_low_pass = ifftshift(Xb_low_pass)
x_low_pass = ifft(X_low_pass)
xb_low_pass = fftshift(x_low_pass.real)

# Function to plot
def plot_signals():
    fig1, axs = plt.subplots(2, 2, figsize=(15, 10))

    time_range = np.arange(-4, 5)
    
    axs[0, 0].stem(time_range, xb, linefmt='blue', markerfmt='bo', basefmt='r-')
    axs[0, 0].set_title('Original Signal')
    axs[0, 0].set_xlabel('Time (s)')
    axs[0, 0].set_ylabel('Amplitude')

    axs[0, 1].stem(time_range, np.abs(Xb), linefmt='red', markerfmt='ro', basefmt='r-')
    axs[0, 1].set_title('Magnitude Spectrum')
    axs[0, 1].set_xlabel('Frequency (Hz)')
    axs[0, 1].set_ylabel('Magnitude')

    axs[1, 0].stem(time_range, Xb_low_pass, linefmt='green', markerfmt='go', basefmt='r-')
    axs[1, 0].set_title('Low-pass Spectrum')
    axs[1, 0].set_xlabel('Frequency (Hz)')
    axs[1, 0].set_ylabel('Magnitude')

    axs[1, 1].stem(time_range, xb_low_pass, linefmt='purple', markerfmt='mo', basefmt='r-')
    axs[1, 1].set_title('Reconstructed Signal')
    axs[1, 1].set_xlabel('Time (s)')
    axs[1, 1].set_ylabel('Amplitude')

    plt.tight_layout()
    plt.show()

plot_signals()

Matlab

close all; clear all;clc;
xb=[1 2 3 4 5 4 3 2 1] % πραγματικό σήμα με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],xb); ylabel('xb');
x=ifftshift(xb) % το σήμα με τις αρνητικές συνιστώσες στο άνω μέρος
X=fft(x) % FFT
Xb=fftshift(X) % το φάσμα με τη dc συνιστώσα στο κέντρο, πραγματικές
% τιμές με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],Xb);ylabel('Xb');
close all; clear all;clc;
Xb=[0 0 1 1 1 1 1 0 0] % φάσμα βαθυπερατού σήματος με άρτια συμμετρία
figure; subplot (2,1,1); plot([-4:4],Xb); ylabel('Xb');
X=ifftshift(Xb) % το φάσμα με τις αρνητικές συνιστώσες στο άνω μέρος
x=ifft(X) % IFFT
xb=fftshift(x) % πραγματικό σήμα με άρτια συμμετρία όπως αναμένεται
subplot (2,1,2); plot([-4:4],xb); ylabel('xb');

Question 3: Modify the previous example so that you start directly with the definition of the spectrum of the low-pass signal \(X\) as expected by the \(ifft\). Write your answer in the text file lab2_nnnnn.txt.

<<<<<<< HEAD

Μέρος 2: Σχεδιασμός και υλοποίηση φίλτρων

Θα ασχοληθούμε με το Παράδειγμα 1.2 της παραγράφου 1.5 του τεύχους Μαθήματος. Το παράδειγμα αυτό παρουσιάζει δύο εναλλακτικές μεθόδους σχεδιασμού FIR φίλτρων:

α) τη μέθοδο των παραθύρων και

β) τη μέθοδο των ισοϋψών κυματώσεων

τις οποίες εφαρμόζει στην περίπτωση βαθυπερατών φίλτρων.

Ο κώδικας σε Matlab περιέχεται στο βιβλίο του μαθήματος στο παράδειγμα 1.2 (Κώδικας 1.3). Στο παράδειγμα, τα φίλτρα δοκιμάζονται σε ένα πραγματικό σήμα, s, το οποίο είναι αποθηκευμένο στο αρχείο sima.mat (binary αρχείο MATLAB). Πρόκειται για ένα σήμα sonar με φάσμα που εκτείνεται μέχρι περίπου τα 4 KHz και συχνότητα δειγματοληψίας Fs=8192 (είναι και αυτή αποθηκευμένη στο αρχείο sima.mat, μαζί με το σήμα). Παρατηρήστε ότι, στον κώδικα του παραδείγματος, προηγείται η φόρτωση του sima.mat (γραμμή 6), γεγονός που επιτρέπει τη χρήση των μεταβλητών s, Fs στο υπόλοιπο τμήμα του (γραμμές 7-44).

Η μέθοδος των παραθύρων

Η μέθοδος των παραθύρων εφαρμόζεται στις γραμμές 8-38 του κώδικα. Πρώτο βήμα αποτελεί ο ορισμός της απόκρισης συχνότητας H ενός ιδανικού βαθυπερατού φίλτρου με συχνότητα αποκοπής \(Fs/8\) (γραμμή 9), μέσω ενός διανύσματος μήκους \(Ν=Fs\), γεγονός που οδηγεί σε ανάλυση συχνότητας \(fο = Fs /N = 1 Hz\). Το διάνυσμα H αποτελείται από μια αλληλουχία μονάδων και μηδενικών που δημιουργούνται από την κλήση των συναρτήσεων \(ones\) και \(zeros\), αντίστοιχα. Για περισσότερες πληροφορίες σχετικά με αυτές τις συναρτήσεις μπορείτε να συμβουλευτείτε την τεκμηρίωση του MATLAB, πληκτρολογώντας \(doc ‘function-name’\) στο παράθυρο εντολών, όπου \(‘function-name’\) το όνομα της συνάρτησης. Τα \(Fs/8\) πρώτα στοιχεία του \(H\), που αντιστοιχούν στις συχνότητες \([0,Fs/8)\), είναι μονάδες, ακολουθούν \(3Fs/4\) μηδενικά και άλλες \(Fs/8\)\( μονάδες που αντιστοιχούν στη ζώνη συχνοτήτων \)[–Fs/8, 0)$.

Επόμενο βήμα είναι ο υπολογισμός του αντίστροφου διακριτού μετασχηματισμού \(Fourier\) \((IDFT)\), που υπολογίζεται από τη συνάρτηση \(ifft\) του MATLAB (γραμμή 12). Ακολουθεί μια αναδιάταξη του αποτελέσματος του αντίστροφου \(DFT\) (γραμμές 13-14), που ισοδυναμεί με ολίσθηση της κρουστικής απόκρισης του φίλτρου κατά το μισό της μήκος. Στη συνέχεια, η κρουστική απόκριση \(h\) περικόπτεται σε μήκη 32+1, 64+1 και 128+1 δειγμάτων (γραμμές 15-17). Με τη βοήθεια του \(wvtool\) (Window Visualization Tool), απεικονίζονται στο ίδιο διάγραμμα η κρουστική απόκριση και η απόκριση συχνότητας του βαθυπερατού φίλτρου και για τα 3 παραπάνω μήκη. Το \(wvtool\) είναι ένα παρόμοιο εργαλείο με το \(wintool\) και χρησιμεύει για την οπτικοποίηση παραθύρων στο πεδίο του χρόνου και της συχνότητας. Παρατηρήστε ότι όσο μεγαλύτερο το μήκος του φίλτρου τόσο μικρότεροι είναι οι πλευρικοί λοβοί στην απόκριση συχνότητας. Για να μειωθούν ακόμα περισσότερο αυτοί οι πλευρικοί λοβοί και οι επιπτώσεις του ορθογωνικού παραθύρου, κατασκευάζονται παράθυρα \(Hamming\) και \(Kaiser\) (μέσω των συναρτήσεων \(hamming\) και \(kaiser\) αντίστοιχα, γραμμές 24-25) τα οποία εφαρμόζονται στο φίλτρο μήκους 64+1 σημείων \(h64\) (γραμμές 27-30). Παρατηρήστε μέσω του \(wvtool\) (γραμμή 31) τη σαφώς χαμηλότερη στάθμη των πλευρικών λοβών στην απόκριση συχνότητας των παραθύρων \(Hamming\) και \(Kaiser\) σε σχέση με το ορθογωνικό. Τέλος, το σήμα s φιλτράρεται με καθένα από τα τρία φίλτρα (ορθογωνικό, \(Hamming\), \(Kaiser\)), χρησιμοποιώντας τη συνάρτηση \(conv()\) που υπολογίζει τη συνέλιξη μεταξύ του σήματος και της κρουστικής απόκρισης του εκάστοτε φίλτρου (γραμμές 33-38). Παράλληλα, υπολογίζεται και σχεδιάζεται η πυκνότητα φάσματος ισχύος κατά \(Welch\) με τη βοήθεια της συνάρτησης \(pwelch\). Το αποτέλεσμα δείχνει εμφανώς την ραγδαία μείωση της φασματικής ισχύος για συχνότητες άνω της \(Fs/8\).

Μέθοδος ισοϋψών κυματώσεων Η μέθοδος των ισοϋψών κυματώσεων εφαρμόζεται στις γραμμές 40-44 του κώδικα για την κατασκευή ενός φίλτρου \(Parks-McClellan\). Η κατασκευή του φίλτρου πραγματοποιείται με την κλήση της συνάρτησης \(firpm\) (γραμμή 41). Για τον ορισμό των παραμέτρων εισόδου της \(firpm\) συμβουλευτείτε την τεκμηρίωση του MATLAB. Η εφαρμογή του φίλτρου στο σήμα γίνεται όπως και στην προηγούμενη μέθοδο, με χρήση της συνάρτησης conv (γραμμή 43).

=======

Part 2: Design and implementation of filters

We will deal with Example 1.2 of paragraph 1.5 of the course’s issue. This example presents two alternative methods for designing FIR filters: a) the window method, and b) the equal ripple method, which it applies in the case of low-pass filters.

Tip

The code in Matlab is included in the course’s book in Example 1.2 (Code 1.3).

In the example, the filters are tested on a real signal, s, which is stored in the file sima.mat (a binary MATLAB file). This is a sonar signal with a spectrum extending up to about 4 KHz and a sampling frequency of Fs=8192 (it is also stored in the sima.mat file, along with the signal). Note that, in the code of the example, the loading of sima.mat (line 6) precedes, allowing the use of the variables s, Fs in the rest of it (lines 7-44).

The Window Method

The window method is applied in lines 8-38 of the code. The first step is the definition of the frequency response H of an ideal low-pass filter with a cutoff frequency of \(Fs/8\) (line 9), using a vector of length \(N=Fs\), leading to a frequency resolution of \(fo = Fs /N = 1 Hz\). The H vector consists of a sequence of ones and zeros created by calling the \(ones\) and \(zeros\) functions, respectively. For more information about these functions, you can consult MATLAB documentation by typing \(doc\) \(‘function-name’\) in the command window, where \(‘function-name’\) is the name of the function. The first \(Fs/8\) elements of \(H\), corresponding to frequencies \([0,Fs/8)\), are ones, followed by \(3Fs/4\) zeros and another \(Fs/8\) ones corresponding to the frequency band \([–Fs/8, 0)\).

The next step is the calculation of the Inverse Discrete Fourier Transform \(IDFT\), computed by MATLAB’s \(ifft\) function (line 12). This is followed by a rearrangement of the inverse \(DFT\) result (lines 13-14), equivalent to shifting the filter’s impulse response by half its length. Then, the impulse response \(h\) is truncated to lengths of 32+1, 64+1, and 128+1 samples (lines 15-17). With the help of \(wvtool\) (Window Visualization Tool), the impulse response and the frequency response of the low-pass filter for the above lengths are displayed on the same diagram. The \(wvtool\) is a tool similar to \(wintool\) and is used for visualizing windows in both time and frequency domains. Notice that the longer the filter length, the smaller the side lobes in the frequency response. To further reduce these side lobes and the effects of the rectangular window, \(Hamming\) and \(Kaiser\) windows are constructed (via the \(hamming\) and \(kaiser\) functions, respectively, lines 24-25) and applied to the 64+1 point length filter \(h64\) (lines 27-30). Through the \(wvtool\) (line 31), observe the significantly lower level of side lobes in the frequency response of the \(Hamming\) and \(Kaiser\) windows compared to the rectangular one. Finally, the signal s is filtered with each of the three filters (rectangular, \(Hamming\), \(Kaiser\)), using the \(conv()\) function that computes the convolution between the signal and the impulse response of each filter (lines 33-38). Concurrently, the power spectrum density by \(Welch\) is calculated and plotted using the \(pwelch\) function. The result clearly shows a sharp reduction in spectral power for frequencies above \(Fs/8\).

Equal Ripple Method

The equal ripple method is applied in lines 40-44 of the code for the construction of a \(Parks-McClellan\) filter. The filter is built by calling the \(firpm\) function (line 41). For the definition of the input parameters of \(firpm\), consult MATLAB documentation. The filter’s application to the signal is done in the same way as the previous method, using the \(conv\) function (line 43).

>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

Code

Python

import numpy as np
from scipy.signal import remez, convolve, welch
import matplotlib.pyplot as plt
import scipy.io
import sounddevice as sd

#6
with open('../files/sima.txt') as f:
    s = [float(x) for x in f]
s=np.array(s)
Fs=8192


# 7
f, Pxx = scipy.signal.welch(s, Fs)
plt.figure(); plt.semilogy(f, Pxx); plt.grid(True); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]')

# 8-9
N = Fs // 2  # Adjusted for Python indexing, assuming Fs is even
H = np.concatenate((np.ones(N//4), np.zeros(N//2), np.ones(N//4)))

# 10-14
h = np.fft.ifft(H, n=N).real
h = np.fft.fftshift(h)  # Equivalent to MATLAB's circshift for centering the impulse response

# 15-17
h32 = h[N//2-16:N//2+17]
h64 = h[N//2-32:N//2+33]
h128 = h[N//2-64:N//2+65]

# 18 - Stem plot for h64
plt.figure()
plt.stem(range(len(h64)), h64, basefmt=" ")
plt.title('Stem Plot of h64')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.grid()

# 19 - Frequency response of h64
w, h = scipy.signal.freqz(h64)
plt.figure()
plt.plot(w / np.pi * (Fs / 2), np.abs(h))  # Normalizing frequency to Hz
plt.title('Frequency Response of h64')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid()

# 20 - For visualizing the frequency responses of h32, h64, and h128 together
plt.figure(figsize=(14, 4))

for filt, label in zip([h32, h64, h128], ['h32', 'h64', 'h128']):
    w, h = scipy.signal.freqz(filt)
    plt.plot(w / np.pi * (Fs / 2), 20 * np.log10(abs(h)), label=label)

plt.title('Frequency Responses of h32, h64, h128')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.grid()
plt.grid(True)
plt.show()

# 22-26
wh = np.hamming(len(h64))
wk = np.kaiser(len(h64), 5)
plt.figure(); plt.plot(wh, 'b', label='Hamming'); plt.plot(wk, 'r', label='Kaiser'); plt.legend(); plt.grid()

# 27-30
h_hamming = h64 * wh
h_kaiser = h64 * wk

# 31 `wvtool` equivalent visualization
# Plot the impulse responses
plt.figure(figsize=(10, 6))
markers, stemlines, baseline = plt.stem(range(len(h64)), h64, linefmt='C0-', markerfmt='C0o', basefmt="C0", label='h64')
plt.setp(markers, markersize = 5)  # Adjusting marker size for better visibility
plt.setp(stemlines, linestyle = '-')  # Adjusting stem line style

markers, stemlines, baseline = plt.stem(range(len(h_hamming)), h_hamming, linefmt='C1-', markerfmt='C1o', basefmt="C1", label='h_hamming')
plt.setp(markers, markersize = 5)
plt.setp(stemlines, linestyle = '-')

markers, stemlines, baseline = plt.stem(range(len(h_kaiser)), h_kaiser, linefmt='C2-', markerfmt='C2o', basefmt="C2", label='h_kaiser')
plt.setp(markers, markersize = 5)
plt.setp(stemlines, linestyle = '-')

plt.title('Impulse Responses')
plt.xlabel('Samples')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# Plot the frequency responses
plt.figure(figsize=(14, 4))
for filt, label, color in zip([h64, h_hamming, h_kaiser], ['h64', 'h_hamming', 'h_kaiser'], ['C0', 'C1', 'C2']):
    w, h = scipy.signal.freqz(filt)
    plt.plot(w / np.pi * (Fs / 2), 20 * np.log10(abs(h)), label=label, color=color)

plt.title('Frequency Responses')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.grid()
plt.show()

# 32-38 Filtering and plotting PSD
y_rect = np.convolve(s, h64, mode='same')
f, Pxx = scipy.signal.welch(y_rect, Fs)
plt.figure();plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h64')

y_hamm = np.convolve(s, h_hamming, mode='same')
f, Pxx = scipy.signal.welch(y_hamm, Fs)
plt.figure();plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h_hamming')

y_kais = np.convolve(s, h_kaiser, mode='same')
f, Pxx = scipy.signal.welch(y_kais, Fs)
plt.figure(); plt.grid(True); plt.semilogy(f, Pxx); plt.xlabel('Frequency [Hz]'); plt.ylabel('PSD [V**2/Hz]'); plt.title('h_kaiser')

# 40-44  Parks-McClellan low-pass filter
freq = np.array([0, 0.10, 0.15, 0.5]) * Fs  # Scale frequencies correctly
gain = [1, 1, 0, 0]

# Call firwin2 with corrected parameters and using the `fs` parameter
hpm = signal.firwin2(65, freq, gain, fs=Fs)
s_pm = np.convolve(s, hpm, 'same')

# Plotting or analyzing the filtered signals (y_rect, y_hamm, y_kais, s_pm) can be done similarly
# For example, to plot the PSD:
f, Pxx_spec = signal.welch(s_pm, Fs, 'flattop', 1024, scaling='spectrum')
plt.figure()
plt.semilogy(f, np.sqrt(Pxx_spec))
plt.xlabel('frequency [Hz]')
plt.ylabel('Linear spectrum [V RMS]')
plt.title('Power spectrum (welch)')
plt.show()

# 45-46 Playing sound (s and s_pm assumed as original and filtered signal)
sd.play(20 * s, Fs); sd.wait()
# sd.play(20 * s_lp, Fs); sd.wait()  # Uncomment and adjust variable name as needed for the filtered signal

Matlab

1. clear all; close all;
2. % Το αρχείο "sima.mat" περιέχει το σήμα s και τη συχνότητα
3. % δειγματοληψίας Fs. Το φάσμα του σήματος εκτείνεται σχεδόν σε όλη την
4. % περιοχή συχνοτήτων μέχρι 4 KHz. Πάνω από 1 KHz, όμως, είναι θόρυβος
5. % και πρέπει να φιλτραριστεί.
6. load sima;
7. figure; pwelch(s,[],[],[],Fs);
8. % Ορίζεται η ιδανική βαθυπερατή συνάρτηση Η, με συχνότ. αποκοπ. Fs/8
9. H=[ones(1,Fs/8) zeros(1,Fs-Fs/4) ones(1,Fs/8)];
10. % Υπολογίζεται η κρουστική απόκριση με αντίστροφο μετασχ. Fourier
11. % Εναλλακτικά, μπορεί να χρησιμοποιηθεί η αναλυτική σχέση Sa(x)
12. h=ifft(H,'symmetric');
13. middle=length(h)/2;
14. h=[h(middle+1:end) h(1:middle)];
15. h32=h(middle+1-16:middle+17);
16. h64=h(middle+1-32:middle+33);
17. h128=h(middle+1-64:middle+65);
18. % figure; stem([0:length(h64)-1],h64); grid;
19. % figure; freqz(h64,1); % σχεδιάζουμε την απόκριση συχνότητας της h64
20. wvtool(h32,h64,h128); % αποκρίσεις συχνότητας των περικομμένων h
21. % Οι πλευρικοί λοβοί είναι υψηλοί!
22. % Πολλαπλασιάζουμε την περικομμένη κρουστική απόκριση με κατάλληλο
23. % παράθυρο. Χρησιμοποιούμε την h64 και παράθυρα hamming και kaiser
24. wh=hamming(length(h64));
25. wk=kaiser(length(h64),5);
26. figure; plot(0:64,wk,'r',0:64,wh,'b'); grid;
27. h_hamming=h64.*wh';
28. % figure; stem([0:length(h64)-1],h_hamming); grid;
29. % figure; freqz(h_hamming,1);
30. h_kaiser=h64.*wk';
31. wvtool(h64,h_hamming,h_kaiser);
32. % Φιλτράρουμε το σήμα μας με καθένα από τα τρία φίλτρα
33. y_rect=conv(s,h64);
34. figure; pwelch(y_rect,[],[],[],Fs);
35. y_hamm=conv(s,h_hamming);
36. figure; pwelch(y_hamm,[],[],[],Fs);
37. y_kais=conv(s,h_kaiser);
38. figure; pwelch(y_kais,[],[],[],Fs);
39. %
40. % Βαθυπερατό Parks-MacClellan
41. hpm=firpm(64, [0 0.10 0.15 0.5]*2, [1 1 0 0]);
42. % figure; freqz(hpm,1);
43. s_pm=conv(s,hpm);
44. figure; pwelch(s_pm,[],[],[],Fs);
45. sound(20*s); % ακούμε το αρχικό σήμα, s
46. sound(20*s_lp); % ακούμε το φιλτραρισμένο σήμα, s_lp

<<<<<<< HEAD

Πειραματισθείτε

1. Τροποποιείστε τον κώδικα στην γραμμή 14 ώστε να επιτύχετε το ίδιο αποτέλεσμα με τη χρήση μιας εκ των συναρτήσεων ifftshift ή fftshift.

2. Τροποποιείστε τον κώδικα ώστε να χρησιμοποιηθεί βαθυπερατό φίλτρο μήκους 140+1 αντί του 64+1. Σχεδιάστε την απόκριση συχνότητας του φίλτρου με σχεδιασμό παραθύρου, όπως στο παράδειγμα, για δύο περιπτώσεις: ορθογωνικού και Hamming.

Ορθογωνικό παράθυρο (απλή περικοπή της h)

=======

Experiment

1. Modify the code on line 14 to achieve the same result using either the ifftshift or fftshift functions.

2. Modify the code to use a low-pass filter of length 140+1 instead of 64+1. Plot the frequency response of the filter with window design, as in the example, for two cases: rectangular and Hamming.

Rectangular Window (simple truncation of h)

>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

%matplotlib inline
import ipywidgets as widgets
from IPython.display import display, clear_output
import matplotlib.pyplot as plt
import numpy as np

# Assuming H is defined and calculated as before
H = np.hstack((np.ones(int(Fs/8)), np.zeros(int(Fs-Fs/4)), np.ones(int(Fs/8))))
h = np.real(np.fft.ifft(H))
middle = int(len(h)/2)
h = np.hstack((h[middle:], h[:middle]))
h32 = h[middle-16:middle+16]
h64 = h[middle-32:middle+32]
h128 = h[middle-64:middle+64]
h160 = h[middle-80:middle+80]
h256 = h[middle-128:middle+128]

h_variants = {
    'h32': h[middle-16:middle+16],
    'h64': h[middle-32:middle+32],
    'h128': h[middle-64:middle+64],
    'h160': h[middle-80:middle+80],
    'h256': h[middle-128:middle+128],
}

output2 = widgets.Output()

# Function to plot
def plot_stem(h_key):
    with output2:
        clear_output(wait=True)  # Clear the previous plots
        h_data = h_variants[h_key]
        x_values = np.arange(len(h_data))
        plt.close('all')  # Close all existing figures
        fig, ax = plt.subplots()
        markerline, stemlines, baseline = ax.stem(x_values, h_data, '-.')
        plt.setp(baseline, 'color', 'k', 'linewidth', 2)
        plt.title('Stem plot of selected filter')
        plt.xlabel('Index')
        plt.ylabel('Amplitude')
        plt.show()

# Setup the widgets
dropdown = widgets.Dropdown(options=list(h_variants.keys()), value='h32', description='Filter:')

# Function that updates the plot based on dropdown
def update_plot(change):
    plot_stem(change['new'])

# Observe dropdown for changes
dropdown.observe(update_plot, names='value')

# Group the dropdown and output widget in a vertical box layout
vbox = widgets.VBox([dropdown, output2])

# Display the VBox
display(vbox)

# Initial call to display the plot
plot_stem(dropdown.value)  # Use the dropdown's value

<<<<<<< HEAD

=======

>>>>>>> 0f723ea9e05ce2442ad422bb888bda5e4a9df6a4

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import ipywidgets as widgets
from IPython.display import display

# Compute frequency responses
freq32, resp32 = signal.freqz(h32)
freq64, resp64 = signal.freqz(h64)
freq128, resp128 = signal.freqz(h128)
freq160, resp160 = signal.freqz(h160)
freq256, resp256 = signal.freqz(h256)

# Function to apply window and compute frequency responses
def compute_filtered_freq_responses(window_type='Rectangular'):
    freqs, resps = {}, {}
    for filt_size in [32, 64, 128, 160, 256]:
        filt_name = f'h{filt_size}'
        filt = eval(filt_name)
        
        if window_type == 'Hamming':
            filt = filt * np.hamming(filt_size)
        elif window_type == 'Kaiser':
            beta = 14  # Example beta value for Kaiser window, adjust as needed
            filt = filt * np.kaiser(filt_size, beta)
            
        freqs[filt_name], resps[filt_name] = signal.freqz(filt)
    return freqs, resps

# Initial computation with Rectangular (no window)
freqs, resps = compute_filtered_freq_responses()

output3 = widgets.Output()

def update_plot(change=None):
    window_type = window_type_dropdown.value
    freqs, resps = compute_filtered_freq_responses(window_type)
    
    with output3:
        clear_output(wait=True)
        plt.figure(figsize=(10, 5))
        for filt in ['h32', 'h64', 'h128', 'h160', 'h256']:
            if checkboxes[filt].value:
                w, h = freqs[filt], resps[filt]
                plt.plot(0.5 * Fs * w / np.pi, 20 * np.log10(np.abs(h)), label=filt)
        plt.title(f'Frequency Response with {window_type} Window')
        plt.xlabel('Frequency (Hz)')
        plt.ylabel('Magnitude (dB)')
        plt.xscale('linear')
        plt.grid(True)
        plt.legend()
        plt.show()

# Create checkboxes for each filter
checkboxes = {f'h{size}': widgets.Checkbox(value=False, description=f'h{size}') for size in [32, 64, 128, 160, 256]}
for cb in checkboxes.values():
    cb.observe(update_plot, names='value')

# Dropdown for selecting the window type
window_type_dropdown = widgets.Dropdown(
    options=['Rectangular', 'Hamming'],
    value='Rectangular',
    description='Window Type:',
)
window_type_dropdown.observe(update_plot, names='value')

# VBox to hold the checkboxes and dropdown
vbox_checkboxes = widgets.VBox(list(checkboxes.values()) + [window_type_dropdown])

# Display the VBox along with the output area
vbox = widgets.VBox([vbox_checkboxes, output3])
display(vbox)

<<<<<<< HEAD

3. Σχεδιάστε φίλτρο Parks-McClellan μήκους 160+1, με τις ίδιες οριακές συχνότητες, όπως στο παράδειγμα (0.1, 0.15). Σχεδιάστε την απόκριση συχνότητας και αυτού του φίλτρου και σχολιάστε τις διαφορές από το ίδιο φίλτρο μήκους 64+1 (του παραδείγματος 1.2). Γράψτε την απάντησή σας στο αρχείο κειμένου lab2_nnnnn.txt.

=======